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The superheated Meissner state in type-I superconductors is studied both analytically and numer- 
ically within the framework of Ginzburg-Landau theory. Using the method of matched asymptotic 
expansions we have developed a systematic expansion for the solutions of the Ginzburg-Landau 
equations in the limit of small re, and have determined the maximum superheating field for the 
existence of the metastable, superheated Meissner state as an expansion in powers of re 1//2 . Our nu- 
merical solutions of these equations agree quite well with the asymptotic solutions for re < 0.5. The 
■ same asymptotic methods are also used to study the stability of the solutions, as well as a modified 

| version of the Ginzburg-Landau equations which incorporates nonlocal electrodynamics. Finally, we 

compare our numerical results for the superheating field for large-re against recent asymptotic results 
for large-re, and again find a close agreement. Our results demonstrate the efficacy of the method 
of matched asymptotic expansions for dealing with problems in inhomogeneous superconductivity 
involving boundary layers. 
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PACS numbers: 74.40. +k, 74.55. +h, 74.60. Ec 



>. . I. INTRODUCTION 

in 



In equilibrium, the superconducting Meissner phase of a bulk type-I superconducting sample exists below a ther- 
modynamic critical field H c (or below H c \ in a type-II superconductor); above this field the sample reverts to the 
normal phase (or the flux lattice phase in a type-II superconductor). However, because the phase transition in both 
iy~j ' cases is first-order it is possible to superheat the Meissner phase and delay the transition to fields well above H c or 
ON H c \. This superheated, metastable Meissner phase is eventually destroyed at a maximum superheating field H^. 

Understanding the origin and stability of the superheated state is the first step in providing a complete description 
of the time-dependent collapse of the Meissner phase, which is important for many applications of type-I supercon- 
ductors. For instance, there have been recent proposals to use superheated type-I superconductors as detectors for 
elementary particles, so that the sample acts as a superconducting "bubble chamber" The passage of a sufficiently 
energetic particle through the sample would initiate the transition to the normal state. Measuring the superheating 
field also provides one of the few methods of experimentally determining the Ginzburg-Landau parameter re in type-I 
Q . superconductors 0. 

The precise value of the maximum superheating field may depend upon extrinsic factors such as defects in the 
• i-H , sample and sample preparation and geometry. If these effects can be minimized then the limit of superheating is 
determined by the boundaries and geometry of the sample. The simplest and most widely studied geometry is a 
superconducting half-space with a magnetic field applied parallel to the surface of the superconductor, which is the 
geometry considered in the remainder of this paper. To model the superconductor we will use the Ginzburg-Landau 
(GL) equations, which provide an accurate description of the surface behavior provided the coherence length £ is large 
compared to microscopic length scales (|] . Previous studies of superheating in type-I superconductors have used a 
variety of heuristic methods to determine the behavior of the GL equations near the surface. Ginzburg Q inferred 
the leading re -1 / 2 dependence of the superheating field from the form of the Ginzburg-Landau equations. Apparently 
unaware of Ginzburg's work, the Orsay group |q] used a variational argument to show that H^/Hc w 2~ 1 / 4 re~ 1 / 2 . 
By combining an ingenious guess for the behavior of the superconducting order parameter near the surface with a 
variational calculation, Parr j^] was able to calculate the next order correction to the Orsay group's result. In addition 
to this analytical work there has also been a great deal of numerical work on solving the GL equations in small-re limit 
|0-1O|, which is reviewed in Ref. [[nj. The numerical results appear to confirm at least some of the analytical work, 
although admittedly over a somewhat restricted range of re. One deficiency common to all of the previous analytical 
approaches is an ad-hoc construction of approximate solutions of the GL equations, leaving us without a procedure 
for systematically improving upon these approximations. In addition, the issue of the stability of the solutions in the 
small- re limit seems not to have been addressed rigorously (for one attempt see |]l2|). 

In this paper we re-examine the problem of superheating in type-I superconductors by using the method of matched 
asymptotic expansions [ p^|Jl~^ 1 to solve the GL equations in the small-re limit. This method was originally developed 
to treat boundary layer problems in fluid mechanics um in a controlled and systematic fashion, and is particularly 
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well suited to the superheating problem, as all of the technical difficulties arise due to a "boundary layer" at the 
surface. Using this method we can calculate the superheating field in the small-re limit as an asymptotic expansion 
in powers of re 1 / 2 , construct uniform asymptotic expansions (i.e., expansions valid for all x as re — > 0) for the order 
parameter and magnetic field, determine the stability of the solutions, and treat nonlocal electrodynamic effects. 
Where appropriate we compare our asymptotic results against numerical solutions of the GL equations, and we 
generally find excellent agreement. We have chosen to present our results in detail, as the methods are probably 
unfamiliar to most physicists. Matched asymptotic expansions have recently been used by Chapman [l5| to study 
superheating in the large-re limit, and our results complement his work. This paper, taken together with Chapman's 
work, demonstrates that these perturbation methods can provide a powerful calculational tool for solving problems 
in inhomogeneous or nonequilibrium superconductivity. 

The remainder of this paper is organized as follows. In Sec. II we describe our numerical methods for solving the GL 
equations. In Sec. Ill wc develop the method of matched asymptotic expansions for the solution of the GL equations 
in the small-re limit, and determine the first three terms in the expansion of fish in powers of re 1 / 2 . In addition, wc 
construct uniform expansions for the order parameter and magnetic field and compare them against our asymptotic 
results. In Sec. IV we examine the second variation of the GL free energy, 8 2 !F, in order to determine the stability 
of our solutions. This is done for both one and two-dimensional perturbations. In Sec. V the method is generalized 
to treat nonlocal electrodynamics. Section VII compares Chapman's asymptotic expansion for H s h for large-re with 
our numerical results, and we find remarkably good agreement. Finally, Sec. VIII is a summary and discussion of our 
results. 

II. NUMERICAL METHODS 



The GL free energy of a superconducting sample occupying the half-space x > is 



F\f,q] = [ d 3 r 

Jx>0 



J_,_.,o I 

re 2 v ' ' ' '2 



(V/) 2 + -(1 - ff + /V + (H a - V x q) 



(2.1) 



where re is the GL parameter, / is the amplitude of the superconducting order parameter, q is the gauge-invariant 
vector potential (h = V x q), and H a is the applied magnetic field. The lengths are in units of the penetration 
depth A and fields are in units of \/2H c . Minimizing this expression with respect to both / and q results in the GL 
equations. In one dimension, with / = f(x) and q = (0, q(x), 0), these equations are 

\f" - Q 2 .f + .f-f = 0, (2.2) 

q" - fq = 0, (2.3) 

h = q 1 . (2.4) 

The task at hand is to solve these equations numerically for a superconducting half-space and to find the largest 
possible applied field (-ff s h) which permits a superconducting solution. To insure that no current passes through the 
boundary at x = and that the sample is totally superconducting infinitely far from the surface, we impose the 
boundary conditions 

/'(0) = 0, asuoo. (2.5) 

Since the field at the surface must equal the applied field H a , and the field infinitely far from the surface must equal 
0, we impose the boundary conditions 

h(0)=H a , q{x)^0 as s -> oo. (2.6) 

For re — > 0, we rescale the equations as x' = rex making the new unit of length the correlation length £. Since £ 3> A 
in this limit, a numerical solution over a domain much larger than £ would insure that the regions of rapid change 
for / and h would be included. (For small re, we find that solving for x' < 500 is sufficient.) In the large re limit, we 
use the rescaled equations again, but we increase the size of the domain depending on the value of re. (The equations 
must be solved for domains as large as x' < 10 4 for values of re ~ I0 3 .) 
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The equations can be solved using the relaxation method Jlq ] . By replacing these ordinary differential equations with 
finite difference equations, one can start with a guess to the solution and iterate using a multi-dimensional Newton's 
method until it relaxes to the true solution. In order to more accurately pick up the detail near the boundary, we 
choose a grid of discrete points with a higher density near x — 0. In particular we choose a density which roughly 
varies as the inverse of the distance from the boundary. (For low k our density, in units of mesh points per coherence 
length, varies approximately from 10 near the boundary to 10 3 at the farthest point from the boundary, while for 
high k it varies from 10 5 to 10~ 2 .) 

ii/sh can be found in the following way. For a given value of k an initial guess is made where there is no applied 
field and the sample is completely superconducting (/ = 1, q = 0, h = 0). The field H a is then "turned up" in small 
increments. For each value of H a a solution is sought using the result from the previous lower field solution as an 
initial guess. Eventually a maximum value for H a is reached, above which one of two things happens: our algorithm 
fails to converge to a solution or it converges to the normal (nonsuperconducting solution). This maximum value of 
H a is the numerical result for H^. Using this algorithm, H s h(K) can be found for a wide range of k's. Each run (for 
a given tt) takes about 60 cpu minutes on an IBM RS 6000/370. We find it sufficient for the purposes of this paper 
to deal with superheating field values for 10~ 3 < n < 10 3 . 

III. ASYMPTOTIC EXPANSIONS FOR SMALL-k 

In this section we will develop an asymptotic expansion for the superheating field for small-K, using the method 
of matched asymptotic expansions ||l3|Jl4}| . For small-K the dominant length scale is the coherence length £, so it is 
natural to have £ serve as our unit of length. This is achieved by rescaling x by k, introducing a new dimensionless 
coordinate x' — kx. The resulting GL equations in these "outer variables" are 

/" -Q 2 .f + f-f = 0, (3.1) 
"V - fq = 0, (3.2) 

h = nq' , (3-3) 

with the primes now denoting differentiations with respect to x' . 

Outer solution. In order to obtain the outer solutions expand /, q, and h in powers of k: 

f = fo + «/i + « 2 / 2 + ■ ■ • , (3.4) 
q = q + nqi + K 2 q 2 + . . . , (3.5) 

h = ho + nh\ + K 2 /i2 + . . . . (3-6) 
Substituting into Eqs. (3.1)-(3.3), at 0(1) we have 

fo - <7o/o + fo ~ fo = 0, (3.7) 

- foqo = 0. (3.8) 



Since we want / — > 1 as x' — > oo, the only possible solution to Eq. (3.8) is qo = 0. We can then immediately integrate 



Eq. (3.7), 



f (x') = tanh (^P- ) , (3.9) 



V2 

with xo an arbitrary constant. To O(k), the outer equations are 

fl - Zqofoqi - qlh + fl - 3/ 2 /l = 0, (3.10) 

- /oV - 2/ogo/i = 0, (3.11) 
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fto = 0. 



(3.12) 



Once again, the only solution to Eq. (3.11) is qi = 0; substituting this into Eq. (3.10), we find /i = Ci/q, with C\ a 
constant: 



= 7f edl 



(3.13) 



(n) 

We can continue in this manner; at every order q n = Q, h n = 0, and /„ = C„/q , with the C„'s constants which are 
determined by matching onto the inner solution. 

Inner solution. The outer solution breaks down within a boundary layer of O(n) near the surface. This suggests 
introducing a rescaled inner coordinate X — x' /k, so that X — 0(1) within the boundary layer. It is also possible 
to rescale / and q, with the hope that this will lead to a tractable inner problem. Such a rescaling must lead to a 
successful matching of the inner and outer solutions; i.e., the inner solutions as X — > oo must match onto the outer 
solutions as x' — > 0. Since /o(0) = tanh(xo / V2) > then assuming that xo ^ we have /o(0) = 0(1), indicating that 
the order parameter should not be rescaled in the inner region; therefore we set f(x') — F(X) in the inner region. 
However, from the outer solution for the vector potential we see that the only constraint on q(X) in the inner region 
is that q(X) — > as X — > oo (presumably exponentially). Therefore, we are free to rescale q by k in the inner region, 
hopefully in a way whic h simplifi es the inner equations. One possibility is q(x') — n~ a Q(X); substituting this into 
the GL equations, Eqs. (3.1)-( |3~3| ), we see that unless 2a is an integer, fractional powers of n will be introduced into 
the inner equations, contradicting our expansion of / and q in integer powers of k in the outer region. Therefore, the 
simplest assumption is that a = 1/2, leading to the following choice for the inner variables: 



x' = kX, f(x') = FIX), q(x') = «T 1/2 QpO, h(x') = H(X) 



In these variables Eqs. (3.1)-(3.3) become 

F" - kQ 2 F + k 2 (F - F 3 ) = 0, 



(3.14) 



(3.15) 



Q" - F 2 Q = 0, 



(3.16) 



k 1 I 2 H = Q', 

where now the primes denote differentiation with respect to X . The boundary conditions are 

F'(0) = 0, H(0) = H a . 
The next step is to expand the inner solutions in powers of k: 

F = F + kF 1 + k 2 F 2 + . . . , 

Q = Qo + kQi + k 2 Q 2 + 



(3.17) 

(3.18) 

(3.19) 
(3.20) 



H = k~ 1/2 H + k 1/2 £/i 



(3.21) 



Note that there is no term of O(l) in the expansion for iJ, since we would be unable to match such a term to the 
outer solution. Using the boundary condition H(0) = H a leads to 



H a = K- 1 l 2 H a {Q) + K 1 ' 2 H 1 {0) + 



Substituting these expansions into Eqs. (3.15)-(3.17), at O(l) we obtain 

Fq = 0, Qo - FqQq = 0, H = Q' . 



(3.22) 



(3.23) 



Solving these equations subject to the boundary conditions ( 3.18 ) (we also need Qo — > as x — * oo in order to match 
onto the outer solution), we obtain 



F {X) = A , Qo{X) = B e 



-AnX 



F (0) = -A B 0l 



(3.24) 
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with Aq and Bq constants. In what follows we will assign F n (0) = A n and Q„(0) = B n for notational simplicity. The 
O(k) equations are 

F[' = Q 2 F , Q'(-F 2 Qx = 2F Q Fx, H 1 = Q' 1 . (3.25) 
Solving with the boundary condition F[(§) = 0, we obtain 

Fx (X) = Ax + [2A X + e- 2A « x - l] , (3.26) 



Qi(X) = e 



-A„X 



Bl - 



\QA 2 



1-e 



-2A X 



AAlX 2 



(3.27) 



Finally, to 0(k 2 ) we have for F 2 



1 B 3 

- A Bx - AxBq. 



Fo 3 



2QoQiFq 



QqFx, 



the solution of which (with -F^O) = 0) is 



F (X) 17 B ° 



1 B 2 A 1 



4 

IBqBx 
2 



lB Bx 
2 



An 



An 



BqBx — 



A 



IB 2 Ax 



32 Al 



3_B4 
32 Al 

8A 2 



X- l -A Q {l-A 2 )X 2 



2 A 



1 B'k n 

X -^x 2 

8A 



3 Bp -aaqX 
128 



(3.28) 



(3.29) 



(3.30) 



The expression for Q 2 is even more unwieldy, and is not needed in what follows. 

Matching. To determine the various integration constants which have been introduced we must match the inner 
solution to the outer solution. Since the outer solution for q is simply q = 0, and all of our inner solutions decay 
exponentially for large X, the matching is automatically satisfied for q, as well as for the magnetic field h. To match 
the inner and outer solutions for the order parameter, we are guided by the van Dyke matching principle fl4|| , which 
states that the m term inner expansion of the n term outer solution should match onto the n term outer expansion 
of the m term inner solution. In our case we will take m = 3 and n = 2. Therefore, write the two term outer solution 
fo(x') + Kfi(x') in terms of the inner variable X, and expand for small K, keeping the first three terms in the expansion 
in powers of k: 



Jo(kX) + nfxinX) ~ tanh (^j + Kscch 2 ^-^0 -L [C x + X] 



-K 2 sech 2 ( — % ) tanh ( — 
\V2j \V2 



(331) 



Next, write the three term inner solution Fq(X) + kFx(X) + k 2 F2(X) in terms of the outer variable x' , and expand 
for small k, this time keeping the first two terms of the expansion: 



W«) + kF^x'/k) + k 2 F,(x'Ik) ~A„ + tS-x 1 -\a„(1- Al)x 12 



(3.32) 



By writing both expressions in terms of x' , and equating the various coefficients of x' and k, we see that the expansions 
do indeed match if we choose 



5 



Aq = tanh 



V2 



(3.33) 



Bo 



_ 2 V4 sech ( ^0 

>V2 



-2 X / 4 (1 - A 2 ) 1 ' 2 , 



(3.34) 



^-4 

4A 



sech^ -i = 



V2J y/2 



An 



(3.35) 



Eliminating Ci 



= 3 £g V2A (l-A 2 Q ) Ci 
1 32 A 2 B V2' 



V2A A! 3B 3 1 1 - j4q 



32 A 2 



2 B n 



Substituting into our expressions for H o (0) and H x (0) from Eqs. ( [Tgg ) and ( fT28| ), we obtain 

#o(0) = 2 1//4 A (1 - Aq) 1 / 2 , 



(3.36) 



(3.37) 



(3.38) 



ffi(0) = 



2 3 / 4 (2A 2 + 14)(1 - Al) 1 ' 2 2 1 / i (2A 2 - 1) 



G4 



A 



(\-A 2 yn 



Al 



(3.39) 



In order to calculate the superheating field (or, more correctly, the maximum superheating field), we need to maximize 
Hq(0) and -Hi(O) with respect to Aq and A±. Maximizing -ffn(O) with respect to Aq, we find that the maximum occurs 
at A* Q = I/a/2, B^ = -2" 1 / 4 , so that H Q (0) = 2~ 3 / 4 . Substituting this result into Hi(0), we find the surprising result 
that the coefficient of A\ is zero, and i?i(0) = 2 3 / 4 15/64. Our superheating field is then 



-3/4^-1/2 



15V2 

/ 

32 



0(k 2 ) 



(3.40) 



In order to determine A\ we need to proceed to a higher order calculation. The method is the same as before, 
although the algebra quickly becomes tedious; we have used the computer algebra system Maple V to organize the 
expansion. The results from a six term inner expansion are summarized in Table |. Including the next order term in 
the expansion in the superheating field, we have 



2 -3/4 K -l/2 



1 H K - 

32 



325 
1024 



K 2 



0(k 3 ) 



(3.41) 



The first term is exactly the result obtained by the Orsay group pLpj , who used a variational argument to obtain their 
result. The second term is identical to the result obtained by Parr |6j. Parr combined an inspired guess for the behavior 
of the order parameter near the surface with a variational calculation in order to obtain his result. It is interesting 
to note that our result for Ax also agrees with Parr's result. The advantage of the method of matched asymptotic 
expansions is that we can make t his ex pansion systematic, and therefore in principle carry out this expansion as far as 
we wish. The third term in Eq. ( 3.41 ) is one of the new results of this paper; the fourth and fifth terms are included 
in Table Q. With the five-term expansion for i? s h it is p ossible to employ resummation techniques to improve the 
expansion. For instance, the [2, 2] Pade approximant |L3|] is 



rrPadc o 



-3/4„-l/2 



1 + 5.4447812 k + 4.2181012 n 2 
1 + 4.7818686 k + 1.3655230 k 2 



(3.42) 



In Fig. [l] we compare the numerically calculated superheating field against the one, two, and three term asymptotic 
expansions. The one term (i.e., the Orsay group) result never seems particularly accurate. There is a marked improve- 
ment with the two term expansion, with the three term expansion offering only a modest additional improvement. 
The [2,2] Pade approximant agrees with the numerical data to within about 1% all the way to K = 1. 
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Uniform solutions. From the inner and outer expansions it is possible to construct uniform solutions, which are 
asymptotically correct for all x as k — > 0. To do this we simply add the inner and outer solutions of a given order, 
which guarantees the correct behavior in the outer region as well as in the boundary layer. However, this would 
produce a result which was 2/ matc h in the matching region, so we need to subtract / ma tch in order to obtain the 
correct behavior in this region. As an example, we will construct the 2-term uniform solution for the order parameter. 
Adding the two-term outer solution, fo(x') + Kfi(x'), to the two-term inner solution, Fo(X) + kFi(X), subtracting 
the solution in the matching region, which is 1/V2+ (v2/4)kX — (15/32)k, and writing the entire combination in 
terms of the original variable x (which is the same as X), we obtain 

/umf,2(z) = tanh f " I - — K sech 2 f ^_ ° j + -e V2x . (3.43) 

As x — > oo, /unif,2(a;) — > 1; also, / un if,2(0) = l/y/2— (7/32)k, as we expect. However, f^ ni{ 2 (0) = (15/64)k 2 , so that 
the zero-derivative boundary condition is only satisfied to 0{n). 

In Fig. | and Fig. | we compare the numerically calculated order parameter and magnetic field with the two term 
outer solutions and the three term inner solutions. The agreement is quite good for k — 0.1, with deviations appearing 
at k = 0.5. These figures also illustrate the existence of a matching region where the inner and outer solutions overlap; 
this region grows as K — * 0. Lastly, we show in Fig. ||how the two term uniform expansion constructed earlier supplies 
a uniform approximation to the order parameter and magnetic field over the whole region. 



IV. STABILITY ANALYSIS OF THE SOLUTIONS 



Having obtained an asymptotic expansion for the superheating field in powers of k 1 / 2 , we now examine the 
stability of the solution with respect to infinitesimal perturbations by studying the second variation of the free energy, 
8 2 T '. Perturbations with 8 2 T > correspond to stable solutions, while those with 5 2 J- < correspond to unstable 
solutions. We will again use the method of matched asymptotics to solve for the eigenfunctions of the linear stability 
operator. We first determine the stability in the simpler one-dimensional situation then we discuss the two-dimensional 
case. 



A. Stability with respect to one-dimensional perturbations 



If we perturb the extremal solution (/, q) of the GL equations by allowing f(x) — > f(x)+f(x) and q(x) — > q(x)+q(x), 
then the second variation of the free energy functional is 



5 2 F = 



dx 



1 X, 



■^f' 2 + (3/ 2 + q Z - l)f 2 + ifqfq + ft + t 



The boundary conditions on / and q should be chosen so as to not perturb / and h at the surface, so that 

/'(0) = q'(0) = 0, /(oo) = g(oo) = 0. 
We can then integrate Eq. fl4.1|) by parts to obtain 



8 2 T 



dx 



k 2 dx 2 

This quadratic form can be conveniently written as 

5 2 T- 

where L\ is the self-adjoint linear operator 



/ ( -4ir + T - :-!/-' - 1 )/ + '/( -^2 + /"' S '/ - J '//'}/ 



dx (/, <z)Li 



1 JL _L ^,2 

dx 



q 2 + 3f 2 - 1 2fq 



2fq 



dx- 



+ f 2 



(4.1) 



(4.2) 



(4.3) 



(4.4) 



(4.5) 
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In order to analyze the stability, expand / and q as 

/ 



(4.6) 



where the c ra 's are real constants, and (/„, q n ) is a normalized eigcnfunction of L\ with eigenvalue E n : 



Li ( l n j = ( ~™ ) • (4.7) 



<7n / V 9 

Then 

<5 2 .F = ^£;„ C 2 . (4.8) 

n 

The second variation iJ 2 ^ 7 ceases to be positive-definite when the lowest eigenvalue first becomes negative, indicating 
that the corresponding solutions (/, q) of the GL equations are unstable. Therefore the entire issue of the stability 
of the solutions has been reduced to finding the eigenvalue spectrum of the linear stability operator L\, which in the 
k — ► limit can be studied using matched asymptotic expansions. 

Outer solution. The outer equations for (/, q) are rescaled with before to yield (we will drop the subscript 

ri for notational convenience) 

- /" + (3/ 2 + q 2 - 1)/ + Ifqq = Ef, (4.9) 



«V + fq + 2fqf = Eq. 



(4-10) 



Expanding /, q, and E in powers of k, and recalling that q = to all orders in k in the outer region, we have at 
leading order 



- fo + (3/o - l)/o = E fo, 



(4.H) 



where fo = tanh y J ■ By changing variables to y = tanh ( - ^ Q j we see that the solution of Eq. ( 4.11 ) is the 
associated Legendre function of the first kind: 



fo(x') = CoPg 



tanh 



x' + x 



(4-12) 



where fi = —^2(2 — Eq) and Co is a constant. The leading order solut ion f or q is qa = 0. 

Inner solution. To obtain the inner equations, we rescale as in Eq. (3.14), with the perturbations rescaled as 

f(x') = F(X), q(x') = k-WQW, (4.13) 

such that 

1 - 1 _ 1 

(4-14) 

F' Z Q + 2FQF = EQ. (4.15) 
To leading order, Fq = 0, so that Fo = aoj with ao a constant. The leading order equation for the variation in Q is 

- Q'o' + 2F Q F + (Fq - E )Q = 0. (4.16) 
The solution which satisfies the boundary condition Q (0) = is 



\p" + (3F 2 + -Q 2 - 1)F + -2FQQ = EF, 



K K 



MX) 



2a A B 
En 



~-AaX 



Ao 



-y/A 2 -E X 



(4.17) 



At O(k) we find 
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(4.18) 



with the solution 

F 1 (X) = a 1 +a B 2 

+a B 2 



E q + AA 2 ^ 2AqX _ 4A% 



-(Ao+^g-EoJX 



So (Ao + y/A* - Eo) V^o - Eo 
4A 3 /E Q 



Eg + 4A 2 

2A E a (A + y/A 2 - E ) y/Al - Eo 



(4.19) 



X. 



We now have enough terms in the inner and outer region for a nontrivial match. 

Matching. We complete the matching of the inner and outer perturbations to obtain the eigenvalue, Eq. Performing 
a two term inner expansion of the one term outer solution, we have 



fo(nX) ~ c 



^sech 2 (x /V2) dP fi Ao) K X 
dA 



(4.20) 



where we have used tanh(xo/v / 2) = Ao- Next, the one term outer expansion of the two term inner solution is 



Fq{x'/k) + kF x {x'/k) ~ a o + 



ao2^(l-^g) 
Eo 



Eo + 4A 2 
2An 



4AI 



(A + ^A 2 -E Q )^A 2 -E 



x', 



(4.21) 



where we have used B = -2 X / 4 (1 - A 2 ,) 1 / 2 . Matching the two expansions using the van Dyke matching principle 
yields 



Co 



«() 



(4.22) 



dP^(Ao) 



P^(Ao) dAo 



2 

Eo~ 



Ep + 4A 2 
2A 



AAl 



{Ao + ^Al-Eo)^Al-Eo 



(4.23) 



The last equation is a rather complicated implicit equation for the eigenvalue Eo(Aq), which generally must be solved 
numerically. However, when Aq = l/y/2 we find Eq = 0, corresponding to the critical case, with E > for Ao > l/y/2. 
The numerical evaluation of Eq. (4.23) is shown in Fig. |^. Therefore, we see that our maximum superheating field 
(at lowest order) corresponds to the limit of metastability for these one-dimensional perturbations. In Fig. |^ we show 
Ao as a function of the lowest order magnetic field at the surface, Ho, from Eq. (3. 35). The stability analysis of this 
section shows that only the upper branch of this double valued function corresponds to solutions which are locally 
stable, with the field at the "nose" being the superheating field. 



B. Stability with respect to two-dimensional perturbations 



We next turn to the stability of the solutions with respect to two dimensional perturbations. If we perturb the 
extremal solution (/, q) of the GL equations by allowing / — ► / + Sf and q — > q + Sq, then the second variation of 
the free energy functional is 



5 2 F 



dx dy 



:(WSf) 2 + 4/(<J/)q • <5q + / 2 (<5q) 2 + (3/ 2 + q 2 - 1)(<5/) 2 + (V X <5q) 2 



(we neglect perturbations along the z-direction) . Expanding in Fourier modes with respect to y JsJ , 
Sf(x,y) = f(x)cosky, 6q x {x, y) = q x (x) sin ky, 5q y (x,y) = q y {x) cosfcy, 



(4.24) 



(4.25) 



substituting into Eq. (4.24), recalling that q = (0, q{x), 0), and integrating over y, we obtain (up to a multiplicative 
constant) 



b 2 T 



dx 



\f' 2 + (3/ 2 + q 2 + \k 2 - l)/ 2 + Afqfq y + f 2 (f x + q 2 ) + (q' - kq x ) 2 



(4.26) 







By integrating by parts and using the boundary conditions, Eq. (|4.2|), we can cast this functional into the form 



5 2 T = 



dx(f,q y ,q x )L 2 




(4.27) 



where the self-adjoint linear operator L 2 is given by 




1 d z 

dx' 2 







■A + 

d 

^ dx 







dx J 


% 


p + ej 





(4.28) 



As in the previous section, we want to determine the eigenvalue spectrum of this operator. We are primarily interested 
in the effects of long-wavelength perturbations (i.e., k — > 0), so we rescale k as k = nk' . Then the eigenvalue equations 
in terms of the outer coordinate x' = kx are (dropping the prime on k from now on) 



/" + (3/ 2 + q 2 - 1 + k 2 )f + 2fqq = Ef, 



* 2 q'y + f\ + 2/<?/~ - K 2 kq' x = Eq y , 



K 2 kq' y + (.T + K Z k l )q x - Eq x 



By using the last equation we may eliminate q x from Eq. ( 4.30 ), which becomes 

/ 2 



dx 



E 



f 2 + K 2 k 2 



E 



+ f 2 q y + 2}qf = Eq y 



(4.29) 
(4.30) 
(4.31) 

(4.32) 



For k — Eqs. ( 4.29 ) and ( 1.32 ) reduce to the one-dimensional perturbation equations of the last section, Eqs. (4.E) 
and (4.1C); for E — they reduce to the Euler-Lagrange equations derived by Kramer [||. 

The perturbation equations (4.29) and ( 4.3 2| ) may be solved by the method of matched asymptotic expansions, just 
as in the one-dimensional case. The derivation of the eigenvalue condition is essentially identical, with the final result 
that 



1 dPg(A ) _ 2 



P^Ao) dA 



E a 



AA 2 



AAl 



2An 



(A + ^A 2 -Eo)^A 2 -E Q 



(4.33) 



where now /i — — \/2(2 + Eq — k 2 ). The eigenvalue Eq (k) is plotted in Fig. for several different values of Aq . For 
A® > l/v2j Eo(k) > for all k, while for Aq < l/v2 there exists a band of long-wavelength perturbations for which 
Eo{k) < 0. In all cases the most unstable modes are at k = 0, i.e., the one-dimensional perturbations are the least 
stable. This is in contrast to the large-K limit, where the most unstable mode occurs for k ^ 0,iL^5| . 



V. NONLOCAL EFFECTS AS k 







In the previous sections we have studied superheating in type-I superconductors starting from the conventional GL 
equations, which assume a local relationship between the current and the vector potential. However, in very clean 
type-I superconductors nonlocal effects are often i mpo rtant (in the Pippard limit; see Ref. jlll). We can model these 
effects by replacing the second GL equation, Eq. (ft.2|), by a nonlocal equation of the form 



(5.1) 



K 2 q" - K(x- x')f 2 (x')q{x') dx' = 0, 
^0 

where K (x — x') is a kernel whose Fourier transform K{k) behaves as 

K(k) : 



A 2 /A| 
a/\k\ 



(local limit); 
(extreme anomalous limit), 



(5.2) 
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with Al the London penetration depth and a a constant |lq] . For A ~ A;, we recover the local limit considered in the 
previous sections of this paper. It is still possible to calculate the superheating field in this nonlocal limit using the 
method of matched asymptotic expansions. Indeed, the prescription is the same as for the local case discussed above; 
we only need to solve a slightly more complicated inner problem. In this section we will calculate the leading order 
superheating field in the nonlocal limit, in order to further illustrate the power and flexibility of our method. 

Outer solution. The outer solution is the same as before; the vec tor potent ial is zero to all orders, and the first two 
terms in the expansion for the order parameter are given by Eqs. (3.9) a nd ( 3T3). 

Inner solution. In the inner region we rescale the variables as in Eq. (3.14). In terms of these variables Eq. (5.1) 
becomes 



/>oo 

Q" - K{X - X')F 2 (X')Q(X')dX' = 0. 
Jo 



(5.3) 



We need to solve this equation, along with the first GL equation, Eq. (3.15), perturbatively in k. Expanding F and 
Q as in Eqs. (3.19)— (3.21), we obtain F (X) = A , as before, and 

/>OC 

Q'o -Al K(X - X')Q (X') dX' = 0. (5.4) 
Jo 

This is an integral equation of the Wiener-Hopf type |l9[ ]. To solve, we Fourier transform, introducing 

/>OC 

Q+(k)= / Q a (X)e lkX dX. (5.5) 
Jo 

After Fourier transforming the integral equation, we perform a Wiener-Hopf factorization J19|, with the result that 



Q+(k)=B, 



iip(k) 



where Bq — Qq(0) is a constant, and 



m = - 

7T 



hi 



J [k? + AlK{k)]W 
k 2 + A 2 K(x)~ 



x 2 +A 2 K(k) 



k 2 



dx. 



(5.6) 



(5.7) 



The Fourier transform can be inverted once a particular form for K(k) is specified (although this is unnecessary for 
the calculation of the superheating field; see below). The magnetic field at this order is H (Q) = —A B as before. 
Proceeding to the next order, we have 



F{' = AoQ 2 . 

By applying the boundary condition F{(0) = 0, we find the general solution 



F 1 (X)=A 1 +A 



X 



x 



Qo(y)dy 



X 



vQliv) dy 



(5.8) 



(5.9) 



with Ai = Fi(0) another constant. The equation for Qi(X) is a rather messy inhomogeneous Wiener-Hopf integral 
equation. Fortunately, its solution is not needed for the leading order calculation of the superheating field. 

Matching. We now turn to the matching of the inner and outer solutions. The two term inner expansion of the one 
term outer solution is 



fo(nX) ~ tanh 



Xq 



,V2/ V2 

The one term outer expansion of the two term inner solution is 

F (x'/k) + kFx(x'/k) - A + A 



sech 



X. 



V2 



fo(y)dy)x' 



(5.10) 



(5.11) 



By using the van Dyke matching principle we find Aq = tanh(a;o/v2) as before, and 
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r°° i 

Ao J Q Ql{y)dy = -j={l-Al). (5.12) 

We can use Parseval's identity to express the left hand side of Eq. (5.12) in terms of |Q+(fc)| 2 , and then use Eq. ( |5.6| ) 
to finally arrive at 



ary., PWi ft = 35(1 -4*)- (5.i3) 

To calculate the superheating field we use Eq. ( 5.1 3| ) to express Bq as a function of Aq; we then substitute this 
result into Hq(0) = —AqBq and maximize with respect to Aq in order to determine the lowest order superheating 
field. In the local limit, K(k) = A 2 /A|, and we obtain Aq = l/y/2 and Hq(0) — 2~ 3 / 4 (Al/A), which is the same as 
our previous result when A ~ Xl- In the extreme anomalous limit K(k) — a/\k\, with a = (37r/4)(A 2 £/A|£o) m the 
Pippard theory |l8f| , where £ is the zero temperature coherence length. Performing the integral, we find 

B = -3 3 / 4 2- 3 /V/%- 1/6 (l - Alf/\ (5.14) 

so that 

H o (0) = 3 3 / 4 2- 3 /V/ 6 4 /6 (l - Al) 1 ' 2 . (5.15) 



The maximum occurs at Aq — y / 5/TT, so that Ho(0) = 0.721 a 1 / 6 . Therefore, the superheating field is 

H sh = 0.721 a 1/6 K~ 1/2 + 0(k 1/2 ) (extreme anomalous limit). (5.16) 

The same result has been obtained by Smith et al. |^p[] using an approximation for the order parameter along with a 
variational calculation (in the spirit of method used by the Orsay group [^J). The advantage of our method is that 
it can be systematically improved upon. Although we have not checked the stability in the extreme anomalous limit, 
the procedure should be entirely analogous to that of the previous section. 



VI. LARGE-k RESULTS 



So far we have used the method of matched asymptotics to solve the GL equations in the small-K limit. Chapman 
|l5| has recently used the same method to treat the one dimensional GL equations in the high-K limit. His final result 
for the superheating field is 

H sh = -L + Ck- 4 / 3 + 0(k- 6 / 3 ) (6.1) 

where the constant C is determined from the solution of the second Painleve transcendent; a numerical evaluation 
yields C — 0.326 [^l). The first term was originally derived by Ginzburg and the second term with the unusual 
dependence upon k is the new term. As seen in Fig || the asymptotic and numerical results agree very well. It turns 
out, however, that the calculated H s h is not actually the superheating field, since the one dimensional solution in the 
large-K limit is unstable with respect to two-dimensional perturbations p^|,pp^]; these instabilities occur at at the 
smaller field — \fB/3^/2 = 0.527. This situation is quite different from that of the small- n limit in which our 
stability calculation (Section IV) found the limit of stability to be right at H^. 



VII. DISCUSSION 



In this paper, the one-dimensional GL equations are solved analytically and numerically for a semi-infinite super- 
conducting sample in the small-K limit in order to determine the maximum superheating field -ff s h- We have used the 
method of matched asymptotic expansions to construct for the first time a systematic perturbative solution of the 
Ginzburg-Landau equations, the results of which agree quite closely with our numerical solutions. The same method 
has been used to determine the stability of these solutions with respect to both one- and two-dimensional infinitesimal 
fluctuations; our analysis shows that two dimensional fluctuations do not lead to any additional destabilizing effects, 
in contrast to the situation in the large-K limit. With little modification this method can also be adapted to treat 
nonlocal electrodynamic effects. Finally, our numerical results for large-K compare well with Chapman's asymptotic 
analysis of this regime. Taken collectively, our results demonstrate the effectiveness of the method of matched asymp- 
totic expansions for dealing with boundary layer problems in the theory of superconductivity. We hope that others 
will find useful applications of the methods developed in this paper in treating inhomogeneous superconductors. 
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FIG. 1. A comparison of the numerically calculated superheating field (heavy line) with the three term asymptotic 
expansion for small-K, and the [2, 2] Pade approximant. The one-term expansion due to the Orsay group deviates systematically 
from the calculated superheating field. The two- and three-term expansions provide a marked improvement over the one-term 
expansion. 

FIG. 2. A comparison of the three term inner and outer solutions for the order parameter and the magnetic field with the 
numerical solution for « = 0.1. The asymptotic solutions approximate the computed values only in the appropriate regions. 
The matching region where the inner and outer meet is O(k) as can be estimated from the inner solution for /. 

FIG. 3. The same as Fig. | for k = 0.5. 
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FIG. 4. A comparison of the two-term uniform solution for the order parameter, f U nU,2{x) (dashed line), with the numerical 
solution (solid line) at k = 0.5. The disagreement of the uniform solution with the boundary condition at x — is of order k 2 . 



FIG. 5. The stability eigenvalue E(Ao), with Ao the value of the order parameter at the surface at leading order. We see 
that E > for Ao > 1/ \/2, indicating locally stable solutions. 

FIG. 6. The order parameter at the surface, Ao, as a function of the field at the surface, Ho, at leading order. The stability 
analysis shows that only the upper branch corresponds to locally stable solutions. The field at the "nose" is the limit of stability, 
and corresponds to the superheating field H = 2~ 3 / 4 = 0.595. 

FIG. 7. The stability eigenvalue E(k) for two-dimensional perturbations of wavenumber k, for several different values of Ao- 
For Ao > 1/ \/2 the eigenvalue is stable for all wavenumbers, while for Ao < 1/ \/2 there exists a band of wavenumbers for which 
the solution is unstable. 



FIG. 8. The numerically calculated superheating field for large-K (solid line), compared with the two-term asymptotic 
expansion derived by Chapman (dashed line). The slope of the dashed line is —4/3. 

TABLE I. Summary of the results of the small- k expansion for the superheating field. Here A n is the value of the order 
parameter F(X) at the surface at n-th order, B n is the value of the vector potential Q(X) at n-th order, C„ is the coefficient 
of the n-th term in the outer expansion of the order parameter, and H n (0) is the n-th order term in the expansion of the 
superheating field. 



n An En Cn gn(0) 



2- 1 / 2 -2- 1 / 4 1 2- 3 / 4 

1 -7/32 -(9/16)2 1 / 4 -(15/16)2 1 / 2 (15/64)2 3 / 4 

2 -(17/1024)2 1/2 (159/2048)2 3/4 225/256 -(325/2048)2 1/4 

3 3211/16384 -(745/4096)2 1/4 -(1125/4096)2 1/2 (14191/65536)2 3/4 

4 -(623575/1572864)2 1/2 (16223049/20971520)2 3/4 16875/131072 -(78495727/62914560)2 1/4 
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